home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 40 / Amiga Format CD40 (1999-05-11)(Future Publishing)(GB)(Track 1 of 3)[!][issue 1999-06].iso / -readerstuff- / paul_qureshi / source / bres.c next >
C/C++ Source or Header  |  1999-03-27  |  11KB  |  468 lines
  1.  
  2. --------------------------------------------------------------------------------
  3. 3d Bresenham Line Drawing  (Why I don't know)
  4.  
  5. void
  6. bresenham_linie_3D(x1, y1, z1, x2, y2, z2)
  7.     int                 x1, y1, z1, x2, y2, z2;
  8. {
  9.     int                 i, dx, dy, dz, l, m, n, x_inc, y_inc, z_inc, 
  10.             err_1, err_2, dx2, dy2, dz2;
  11.     int                 pixel[3];
  12.  
  13.     pixel[0] = x1;
  14.     pixel[1] = y1;
  15.     pixel[2] = z1;
  16.     dx = x2 - x1;
  17.     dy = y2 - y1;
  18.     dz = z2 - z1;
  19.     x_inc = (dx < 0) ? -1 : 1;
  20.     l = abs(dx);
  21.     y_inc = (dy < 0) ? -1 : 1;
  22.     m = abs(dy);
  23.     z_inc = (dz < 0) ? -1 : 1;
  24.     n = abs(dz);
  25.     dx2 = l << 1;
  26.     dy2 = m << 1;
  27.     dz2 = n << 1;
  28.  
  29.     if ((l >= m) && (l >= n)) {
  30.         err_1 = dy2 - l;
  31.         err_2 = dz2 - l;
  32.         for (i = 0; i < l; i++) {
  33.             PUT_PIXEL(pixel);
  34.             if (err_1 > 0) {
  35.                 pixel[1] += y_inc;
  36.                 err_1 -= dx2;
  37.             }
  38.             if (err_2 > 0) {
  39.                 pixel[2] += z_inc;
  40.                 err_2 -= dx2;
  41.             }
  42.             err_1 += dy2;
  43.             err_2 += dz2;
  44.             pixel[0] += x_inc;
  45.         }
  46.     } else if ((m >= l) && (m >= n)) {
  47.         err_1 = dx2 - m;
  48.         err_2 = dz2 - m;
  49.         for (i = 0; i < m; i++) {
  50.             PUT_PIXEL(pixel);
  51.             if (err_1 > 0) {
  52.                 pixel[0] += x_inc;
  53.                 err_1 -= dy2;
  54.             }
  55.             if (err_2 > 0) {
  56.                 pixel[2] += z_inc;
  57.                 err_2 -= dy2;
  58.             }
  59.             err_1 += dx2;
  60.             err_2 += dz2;
  61.             pixel[1] += y_inc;
  62.         }
  63.     } else {
  64.         err_1 = dy2 - n;
  65.         err_2 = dx2 - n;
  66.         for (i = 0; i < n; i++) {
  67.             PUT_PIXEL(pixel);
  68.             if (err_1 > 0) {
  69.                 pixel[1] += y_inc;
  70.                 err_1 -= dz2;
  71.             }
  72.             if (err_2 > 0) {
  73.                 pixel[0] += x_inc;
  74.                 err_2 -= dz2;
  75.             }
  76.             err_1 += dy2;
  77.             err_2 += dx2;
  78.             pixel[2] += z_inc;
  79.         }
  80.     }
  81.     PUT_PIXEL(pixel);
  82. }
  83.  
  84. -------------------------------------------------------------------------------
  85. Bresenham's algorithm for ellipses
  86.  
  87. void symmetry( int x, int y )
  88. {
  89.     PUT_PIXEL ( -x, -y );  PUT_PIXEL ( +x, -y );
  90.     PUT_PIXEL ( -x, +y );  PUT_PIXEL ( +x, +y );
  91. }
  92.  
  93. void bresenham_ellipse( int a, int b )
  94. {
  95.     int x,y,a2,b2, S, T;
  96.     
  97.     a2 = a*a;
  98.     b2 = b*b;
  99.     x = 0;
  100.     y = b;
  101.     S = a2*(1-2*b) + 2*b2;
  102.     T = b2 - 2*a2*(2*b-1);
  103.     symmetry(x,y);
  104.     do {
  105.         if (S<0) {
  106.             S += 2*b2*(2*x+3);
  107.             T += 4*b2*(x+1);
  108.             x++;
  109.         } else if (T<0) {
  110.             S += 2*b2*(2*x+3) - 4*a2*(y-1);
  111.             T += 4*b2*(x+1) - 2*a2*(2*y-3);
  112.             x++;
  113.             y--;
  114.         } else {
  115.             S -= 4*a2*(y-1);
  116.             T -= 2*a2*(2*y-3);
  117.             y--;
  118.         }
  119.         symmetry(x,y);
  120.     } while (y>0);
  121. }
  122.  
  123. -------------------------------------------------------------------------------
  124. //
  125. // CONIC  2D Bresenham-like conic drawer.
  126. //   CONIC(Sx,Sy, Ex,Ey, A,B,C,D,E,F) draws the conic specified
  127. //   by A x^2 + B x y + C y^2 + D x + E y + F = 0, between the
  128. //   start point (Sx, Sy) and endpoint (Ex,Ey).
  129.  
  130. // Author: Andrew W. Fitzgibbon (andrewfg@ed.ac.uk),
  131. //     Machine Vision Unit,
  132. //         Dept. of Artificial Intelligence,
  133. //     Edinburgh University, 5 Forrest Hill, EH1 2QL, UK
  134. //     
  135. // Date: 31-Mar-94
  136.  
  137. #include <stdlib.h>
  138. #include <stdio.h>
  139. #include <math.h>
  140.  
  141. static int DIAGx[] = {999, 1,  1, -1, -1, -1, -1,  1,  1};
  142. static int DIAGy[] = {999, 1,  1,  1,  1, -1, -1, -1, -1};
  143. static int SIDEx[] = {999, 1,  0,  0, -1, -1,  0,  0,  1};
  144. static int SIDEy[] = {999, 0,  1,  1,  0,  0, -1, -1,  0};
  145. static int BSIGNS[] = {99, 1,  1, -1, -1,  1,  1, -1, -1};
  146.  
  147. int   debugging = 1;
  148.  
  149. struct ConicPlotter {
  150.   virtual void plot(int x, int y);
  151. };
  152.  
  153. struct DebugPlotter : public ConicPlotter {
  154.   int xs;
  155.   int ys;
  156.   int xe;
  157.   int ye;
  158.   int A;
  159.   int B;
  160.   int C;
  161.   int D;
  162.   int E;
  163.   int F;      
  164.  
  165.   int octant;
  166.   int d;
  167.  
  168.   void plot(int x, int y);
  169. };
  170.  
  171. void DebugPlotter::plot(int x, int y)
  172. {
  173.   printf("%3d %3d\n",x,y);
  174.  
  175.   if (debugging) {
  176.     // Translate start point to origin...
  177.     float tF = A*xs*xs + B*xs*ys + C*ys*ys + D*xs + E*ys + F;
  178.     float tD = D + 2 * A * xs + B * ys;
  179.     float tE = E + B * xs + 2 * C * ys;
  180.   
  181.     float tx = x - xs + ((float)DIAGx[octant] + SIDEx[octant])/2;
  182.     float ty = y - ys + ((float)DIAGy[octant] + SIDEy[octant])/2;
  183.     // Calculate F
  184.     
  185.     float td = 4*(A*tx*tx + B*tx*ty + C*ty*ty + tD*tx + tE*ty + tF);
  186.     
  187.     fprintf(stderr,"O%d ", octant);
  188.     if (d<0)
  189.       fprintf(stderr," Inside "); 
  190.     else 
  191.       fprintf(stderr,"Outside "); 
  192.     float err = td - d;
  193.     fprintf(stderr,"Real(%5.1f,%5.1f) = %8.2f Recurred = %8.2f err = %g\n", 
  194.         tx, ty, td/4, d/4.0f, err);
  195.     if (fabs(err) > 1e-14)
  196.       abort();
  197.   }
  198.   
  199. }
  200.  
  201. inline int odd(int n)
  202. {
  203.   return n&1;
  204. }
  205.  
  206. inline int abs(int a)
  207. {
  208.   if (a > 0)
  209.     return a;
  210.   else
  211.     return -a;
  212. }
  213.     
  214. int getoctant(int gx, int gy)
  215. {
  216.   // Use gradient to identify octant.
  217.   int upper = abs(gx)>abs(gy);
  218.   if (gx>=0)                // Right-pointing
  219.     if (gy>=0)              //    Up
  220.       return 4 - upper;
  221.     else                //    Down
  222.       return 1 + upper;
  223.   else                  // Left
  224.     if (gy>0)               //    Up
  225.       return 5 + upper;
  226.     else                //    Down
  227.       return 8 - upper;
  228. }
  229.  
  230. int conic(int xs, int ys, int xe, int ye,
  231.       int A, int B, int C, int D, int E, int F,
  232.       ConicPlotter * plotterdata)
  233. {
  234.   A *= 4;
  235.   B *= 4;
  236.   C *= 4;
  237.   D *= 4;
  238.   E *= 4;
  239.   F *= 4;
  240.   
  241.   // Translate start point to origin...
  242.   F = A*xs*xs + B*xs*ys + C*ys*ys + D*xs + E*ys + F;
  243.   D = D + 2 * A * xs + B * ys;
  244.   E = E + B * xs + 2 * C * ys;
  245.   
  246.   // Work out starting octant
  247.   int octant = getoctant(D,E);
  248.   
  249.   int dxS = SIDEx[octant]; 
  250.   int dyS = SIDEy[octant]; 
  251.   int dxD = DIAGx[octant];
  252.   int dyD = DIAGy[octant];
  253.  
  254.   int bsign = BSIGNS[octant];
  255.   int d,u,v;
  256.   switch (octant) {
  257.   case 1:
  258.     d = A + B/2 + C/4 + D + E/2 + F;
  259.     u = A + B/2 + D;
  260.     v = u + E;
  261.     break;
  262.   case 2:
  263.     d = A/4 + B/2 + C + D/2 + E + F;
  264.     u = B/2 + C + E;
  265.     v = u + D;
  266.     break;
  267.   case 3:
  268.     d = A/4 - B/2 + C - D/2 + E + F;
  269.     u = -B/2 + C + E;
  270.     v = u - D;
  271.     break;
  272.   case 4:
  273.     d = A - B/2 + C/4 - D + E/2 + F;
  274.     u = A - B/2 - D;
  275.     v = u + E;
  276.     break;
  277.   case 5:
  278.     d = A + B/2 + C/4 - D - E/2 + F;
  279.     u = A + B/2 - D;
  280.     v = u - E;
  281.     break;
  282.   case 6:
  283.     d = A/4 + B/2 + C - D/2 - E + F;
  284.     u = B/2 + C - E;
  285.     v = u - D;
  286.     break;
  287.   case 7:
  288.     d = A/4 - B/2 + C + D/2 - E + F;
  289.     u =  -B/2 + C - E;
  290.     v = u + D;
  291.     break;
  292.   case 8:
  293.     d = A - B/2 + C/4 + D - E/2 + F;
  294.     u = A - B/2 + D;
  295.     v = u - E;
  296.     break;
  297.   default:
  298.     fprintf(stderr,"FUNNY OCTANT\n");
  299.     abort();
  300.   }
  301.   
  302.   int k1sign = dyS*dyD;
  303.   int k1 = 2 * (A + k1sign * (C - A));
  304.   int Bsign = dxD*dyD;
  305.   int k2 = k1 + Bsign * B;
  306.   int k3 = 2 * (A + C + Bsign * B);
  307.  
  308.   // Work out gradient at endpoint
  309.   int gxe = xe - xs;
  310.   int gye = ye - ys;
  311.   int gx = 2*A*gxe +   B*gye + D;
  312.   int gy =   B*gxe + 2*C*gye + E;
  313.   
  314.   int octantcount = getoctant(gx,gy) - octant;
  315.   if (octantcount <= 0)
  316.     octantcount = octantcount + 8;
  317.   fprintf(stderr,"octantcount = %d\n", octantcount);
  318.   
  319.   int x = xs;
  320.   int y = ys;
  321.   
  322.   while (octantcount > 0) {
  323.     if (debugging)
  324.       fprintf(stderr,"-- %d -------------------------\n", octant); 
  325.     
  326.     if (odd(octant)) {
  327.       while (2*v <= k2) {
  328.     // Plot this point
  329.     ((DebugPlotter*)plotterdata)->octant = octant;
  330.     ((DebugPlotter*)plotterdata)->d = d;
  331.     plotterdata->plot(x,y);
  332.     
  333.     // Are we inside or outside?
  334.     if (d < 0) {            // Inside
  335.       x = x + dxS;
  336.       y = y + dyS;
  337.       u = u + k1;
  338.       v = v + k2;
  339.       d = d + u;
  340.     }
  341.     else {              // outside
  342.       x = x + dxD;
  343.       y = y + dyD;
  344.       u = u + k2;
  345.       v = v + k3;
  346.       d = d + v;
  347.     }
  348.       }
  349.     
  350.       d = d - u + v/2 - k2/2 + 3*k3/8; 
  351.       // error (^) in Foley and van Dam p 959, "2nd ed, revised 5th printing"
  352.       u = -u + v - k2/2 + k3/2;
  353.       v = v - k2 + k3/2;
  354.       k1 = k1 - 2*k2 + k3;
  355.       k2 = k3 - k2;
  356.       int tmp = dxS; dxS = -dyS; dyS = tmp;
  357.     }
  358.     else {              // Octant is even
  359.       while (2*u < k2) {
  360.     // Plot this point
  361.     ((DebugPlotter*)plotterdata)->octant = octant;
  362.     ((DebugPlotter*)plotterdata)->d = d;
  363.     plotterdata->plot(x,y);
  364.     
  365.     // Are we inside or outside?
  366.     if (d > 0) {            // Outside
  367.       x = x + dxS;
  368.       y = y + dyS;
  369.       u = u + k1;
  370.       v = v + k2;
  371.       d = d + u;
  372.     }
  373.     else {              // Inside
  374.       x = x + dxD;
  375.       y = y + dyD;
  376.       u = u + k2;
  377.       v = v + k3;
  378.       d = d + v;
  379.     }
  380.       }
  381.       int tmpdk = k1 - k2;
  382.       d = d + u - v + tmpdk;
  383.       v = 2*u - v + tmpdk;
  384.       u = u + tmpdk;
  385.       k3 = k3 + 4*tmpdk;
  386.       k2 = k1 + tmpdk;
  387.       
  388.       int tmp = dxD; dxD = -dyD; dyD = tmp;
  389.     }
  390.     
  391.     octant = (octant&7)+1;
  392.     octantcount--;
  393.   }
  394.  
  395.   // Draw final octant until we reach the endpoint
  396.   if (debugging)
  397.     fprintf(stderr,"-- %d (final) -----------------\n", octant); 
  398.     
  399.   if (odd(octant)) {
  400.     while (2*v <= k2 && x != xe && y != ye) {
  401.       // Plot this point
  402.       ((DebugPlotter*)plotterdata)->octant = octant;
  403.       ((DebugPlotter*)plotterdata)->d = d;
  404.       plotterdata->plot(x,y);
  405.       
  406.       // Are we inside or outside?
  407.       if (d < 0) {          // Inside
  408.     x = x + dxS;
  409.     y = y + dyS;
  410.     u = u + k1;
  411.     v = v + k2;
  412.     d = d + u;
  413.       }
  414.       else {                // outside
  415.     x = x + dxD;
  416.     y = y + dyD;
  417.     u = u + k2;
  418.     v = v + k3;
  419.     d = d + v;
  420.       }
  421.     }
  422.   }
  423.   else {                // Octant is even
  424.     while ((2*u < k2) && (x != xe) && (y != ye)) {
  425.       // Plot this point
  426.       ((DebugPlotter*)plotterdata)->octant = octant;
  427.       ((DebugPlotter*)plotterdata)->d = d;
  428.       plotterdata->plot(x,y);
  429.       
  430.       // Are we inside or outside?
  431.       if (d > 0) {          // Outside
  432.     x = x + dxS;
  433.     y = y + dyS;
  434.     u = u + k1;
  435.     v = v + k2;
  436.     d = d + u;
  437.       }
  438.       else {                // Inside
  439.     x = x + dxD;
  440.     y = y + dyD;
  441.     u = u + k2;
  442.     v = v + k3;
  443.     d = d + v;
  444.       }
  445.     }
  446.   }
  447.  
  448.  
  449.  
  450.   return 1;
  451. }
  452.  
  453. main(int argc, char ** argv)
  454. {
  455.   DebugPlotter db;
  456.   db.xs = -7;
  457.   db.ys = -19;
  458.   db.xe = -8;
  459.   db.ye = -8;
  460.   db.A = 1424;
  461.   db.B = -964;
  462.   db.C = 276;
  463.   db.D = 0;
  464.   db.E = 0;
  465.   db.F = -40000;
  466.   conic(db.xs,db.ys,db.xe,db.ye,db.A,db.B,db.C,db.D,db.E,db.F, &db);
  467. }
  468.